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ABSTRACT 

We present a new map-making method for CMB measurements. The method is based 
on the destriping technique, but it also utilizes information about the noise spectrum. 
The low-frequency component of the instrument noise stream is modelled as a super- 
position of a set of simple base functions, whose amplitudes are determined by means 
of maximum- likelihood analysis, involving the covariance matrix of the amplitudes. We 
present simulation results with 1// noise and show a reduction in the residual noise 
with respect to ordinary destriping. This study is related to Planck LFI activities. 

Key words: methods: data analysis - cosmology: cosmic microwave background 



1 INTRODUCTION 

Construction of the CMB map from time-ordered data 
(TOD) is an important part of the data analysis of CMB ex- 
periments. Future space missions like Planck 1 present new 
challenges for the data analysis. The amount of data Planck 
produces is far larger than that of any earlier experiments. 

The destriping technique (Burigana et al. 1997, De- 
labrouille 1998, Maino et al. 1999, 2002, Keihanen et al. 
2004) provides an efficient map-making method for large 
data sets. The method is non-optimal in accuracy but fast 
and stable. Other methods have been developed which aim 
at finding the optimal minimum-variance map for Planck - 
like data (e.g. Natoli et al. 2001, Dore et al. 2001, Borrill 
1999). 

In this paper we present a new map-making method for 
CMB measurements, called MADAM (Map-making through 
Destriping for Anisotropy Measurements). The method is 
built on the destriping technique but unlike ordinary destrip- 
ing, it also utilizes information about the noise spectrum. 
The aim is to improve the accuracy of the output map as 
compared to destriping, while still keeping, at least partly, 
the speed and stability of the destriping method. 

The basic idea of the method is the following. The low- 
frequency component of the instrument noise in the TOD is 
modelled as a superposition of simple base functions, whose 
amplitudes are determined by means of maximum-likelihood 
analysis, involving the covariance matrix of the amplitudes. 
The covariance matrix is computed from the noise spectrum, 
assumed to be known. 

This paper is organized as follows. In Section 2 we go 
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through the maximum-likelihood analysis that forms the ba- 
sis of our map-making method. We describe the map- making 
algorithm in Section 3. In Section 4 we consider the covari- 
ance matrix of component amplitudes. Some technical cal- 
culations related to this section are presented in Appenxix 
A. We give results from numerical simulations in Section 5 
and present our conclusions in Section 6. 



2 MAP-MAKING PROBLEM 

2.1 Maximum-likelihood analysis of the destriping 
problem 

In the following we present a maximum-likelihood analysis 
on which our map-making me thod is ba s ed. T he analysis 
is similar to that presented in iKeihane'nl i2004T) , the main 
difference being, that here we include the covariance of the 
correlated noise component. 

We write the time-ordered data (TOD) stream as 



= Pm + n'. 



(1) 



Here the first term presents the CMB signal and the second 
term presents noise. Vector m presents the pixelized CMB 
map and pointing matrix P spreads it into TOD. 

We divide the noise contribution into a correlated noise 
component and white noise, and model the correlated part 
as a linear combination of some orthogonal base functions, 



Fa 



(2) 



Vector a contains the unknown amplitudes of the base func- 
tions and matrix F spreads them into TOD. Each column of 
matrix F contains the values of the corresponding base func- 
tion along the TOD. Assuming that the white noise compo- 
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nent and the correlated noise component are independent, 
the total noise covariance is given by 

C t = (n'(n') T ) = FC a F T + C„ (3) 

where C„ = (nn T ) is the white noise covariance, C a = 
(aa T ) is the covariance matrix for the component amplitudes 
a, and (a;) denotes the expectation value of quantity x. 

Our aim is to find the maximum-likelihood estimate of 
m and a simultaneously, for given data y. We consider the 
likelihood of the data, which by the chain rule of probability 
theory can be written as 

P(y) = P(y|m, a)P(a|m)P(m). (4) 

Here P(a\b) denotes the conditional probability of a under 
condition b. Now -P(m) is constant, since we are treating 
the underlying CMB sky as deterministic (we associate no 
probability distribution to it). The probability distribution 
of a is independent of the map so that P(a|m) = P(a). We 
assume gaussianity and write 

P(a) = (27rdetC a )- 1/2 exp(-ia T C- 1 a). (5) 

With m and a fixed, the likelihood of the data is given by 
the white noise distribution 

P(y|m,a) = (27rdet C n )" 1/2 exp(-in T C- 1 n) (6) 

where now n = y — Fa — Pm. The white noise covariance 
C n is assumed to be diagonal, but not necessarily uniform. 

Maximizing the likelihood Q is equivalent to mimimiz- 
ing the inverse of its logarithm. We obtain the chi-square 
minimization function 

X 2 = -21nP(y) =-21n(P(y|m,a)P(a)) 
= (y-Fa-Pm) T C- 1 (y-Fa-Pm) 

+a T C„" 1 a + constant (7) 

We want to minimize Q with respect to both a and m. 
Minimization with respect to m gives 

m= (P T C- 1 P)- 1 P T C- 1 (y-Fa). (8) 

Substituting Eq. J5J back into Eq. Q we get the chi-square 
into the form 

X 2 = (y - Fa) T Z T C,T 1 Z(y - Fa) + a^C^a, (9) 
where 

Z = I-P(P T C- 1 P)- 1 P T C" 1 . (10) 

Here I denotes the unit matrix. 

We minimize \ 2 with respect to a, to obtain an estimate 
for the amplitude vector a. The solution is given by 

(F T C^ 1 ZF + C- 1 )a = F T C r ; 1 Zy. (11) 

Here we have used the property Z T C^ 1 Z = C^ Z. 

The MADAM algorithm uses the conjugate gradient 
technique to solve vector a from Hill . Note that the ma- 
trix on the left-hand side is symmetric. An estimate for the 
CMB map can then be computed using Eq. JSJ. 

2.2 Comparison to the minimum- variance solution 

If the underlying CMB map is treated as deterministic, noise 
is Gaussian, and its statistical properties are known, the 
optimal minimum- variance map is given by 



m = (P T C t - 1 P)- 1 P T C t - 1 y (12) 

where Ct is the covariance of the noise TOD. 

In the following we show that if the total noise covari- 
ance is of the form ©, the map estimate given by Eqs. ||HJ 
and (li lt equals the minimum- variance map 1121 . 

We develop the solution 11H into Taylor series as 

a = (C a - C a F T C- 1 ZFC a + . . .)F T C- 1 Zy. (13) 

We write y — Fa out and recollapse the resulting expansion, 
to get 

y-Fa= (I + FC a F T C- 1 Z)- 1 y (14) 

We now use Eq. Q and write FC„F T C^ 1 = CtC" 1 - I. 
Using this and writing Z out we arrive at 

P T C; 1 (y-Fa) (15) 
= P T [C 4 - (C.C- 1 - I)P(P T C- 1 P)- 1 P T ]- 1 y 

= [I - p t (c,t 1 - cr 1 )P(P T c- 1 P)- 1 ]- 1 P T cr 1 y- 

In the last equality we have taken Ct out from the left and 
used the identity A(I + BA)" 1 = (I + AB) _1 A, which is 
easily verified by expanding both sides as Taylor series. The 
MADAM solution for the map (Eq. JSJ) then becomes 

m = [P T C.- 1 P - P T (C- 1 - C t - 1 )P]- 1 P T C t - 1 y (16) 

which readily simplifies into 1121 . 

If the chosen base function set accurately models the 
correlated noise component, the CMB map estimate given 
by MADAM equals the minimum- variance solution. This is 
necessarily true at the limit where the number of base func- 
tions L per ring approaches the number of samples n, since 
the base functions then form a complete orthogonal basis. 
In practice, however, it is not possible to use that many 
base functions, since both the required memory and CPU 
time increase with increasing L so that the method finally 
becomes unfeasible. 



3 IMPLEMENTATION 
3.1 Map-making algorithm 

Equations © and 1111 form the basis of the MADAM map- 
making method. In this Section we consider the implemen- 
tation of the method. 

Our starting point is a Planck -like scanning strategy, 
where the detector scans the CMB sky in circles which fall on 
top of each other on consecutive rotations of the instrument. 
In order to reduce the amount of data, the data from con- 
secutive scan circles is averaged, a process called 'coadding'. 
We call one coadded circle a 'ring'. In the nominal Planck 
scanning strategy, the same circle is scanned 60 times be- 
fore repointing. We denote by M the number of coadded 
circles. On each ring, we model the correlated noise com- 
ponent as a linear combination of some simple orthogonal 
arithmetic functions, such as sine and cosine functions or 
Legendre polynomials. In the simplest case we fit only uni- 
form baselines. 

The MADAM algorithm can easily be generalized to a 
scanning strategy with no coadding, by setting the coadding 
factor M equal to one. We consider both types of scanning 
strategy in the simulation section. If no coadding is done, 
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Table 1. Main parameters and symbols used in this paper and 
values used in simulation. 



symbol 


parameter 


value 


n 


number of samples / ring 


4608 


N 


number of rings 


8640 


M 


coadding factor 


60 


L 


number of base functions 


1-65 


fs 


sampling frequency 


76.8 Hz 


/sp 


spin frequency of the spacecraft 


1/60 Hz 


a 


white noise std 


2700 /iK 


/mill 


minimum freq. (noise spectrum) 


10" 5 Hz 


fkn 


knee frequency (noise spectrum) 


0.1 Hz 


7 


slope of the noise spectrum 


-1.0 


m 


pixelized CMB map 




P 


pointing matrix 




Cn 


white noise covariance 




F 


base function matrix 




a 


base function amplitudes 




a 


reference values for a 




C„ 


covariance of amplitudes a 




C* 


fcth component covariance 




b k 


... and its coefficient 




9k 


fcth characteristic frequency 





one can choose any length of data to represent a ring. The 
concept of ring then loses its connection to the spin period 
and its length becomes a freely chosen parameter. 

The most frequently used parameters and symbols of 
this paper have been collected in Tabled 

The MADAM map- making algorithm consists of the fol- 
lowing steps. 

(i) Choose a set of base functions to model the correlated 
noise component and compute the covariance matrix C a of 
their amplitudes. The computation of the covariance matrix 
is discussed in Section 4. 

(ii) Using the conjugate gradient technique, solve a from 
Eq. 1111 . The tricky part here is the evaluation of the term 
C^a, since matrix C a is very large. For Planck -like data 
its dimension typically varies from thousands to hundreds of 
thousands. Fortunately, the matrix has symmetries, which 
allow us to evaluate C^a in a quite efficient manner, as we 
show in Section f3. 21 

(iii) Solve the CMB map according to Eq. JHJ. This means 
simple binning of the destriped TOD into pixels, weighting 
by the white noise covariance. Here we have used HEALPix 
2 (Gorski et al. 1999, 2004) pixelization. 

3.2 Evaluation of C^a 

Conjugate gradient solution of Eq. 11 U requires that we 
evaluate 

x = C" 1 a (17) 

several times for different a. Here a and x are vectors of 
(NL) elements, where N is the number of rings and L is the 
number of base functions per ring. Matrix C a has dimension 
(NL, NL) and is thus expensive to invert. However, C a has 
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symmetries which allow us to perform the inversion quite 
efficiently. 

We use index notation in the following. Evaluation Eq. 
<17> is equivalent to solving xu from 

an = y^Cq.i-i'u'ffi'i'- (18) 

i'V 

Here i, i' label rings and I, I' label base functions. The matrix 
has the symmetry property C a>u nv = C ati i Wl . 

We assume that the properties of the correlated noise 
component do not change with time. Matrix C a then de- 
pends on indices i, i' only through their difference, being 
thus approximately circulant in indices . The matrix can 
be stored as a table of L 2 N elements. 

A general symmetric matrix equation of moderate size 
can be solved by Cholesky decomposition. Crout's algorithm 
to find the decomposition is given e.g. in Press et al. (1992). 
On the other hand, circulant matrix equations can be solved 
by the Fourier technique. 

We solve equation (1181 using a combined technique, 
where we apply Cholesky decomposition to the indices 1,1' , 
and Fourier technique to indices 

We drop indices i,i for a while and introduce the fol- 
lowing notation. We denote by Giy an (N, N) submatrix of 
C a - Now C can be understood as an (L, L) matrix whose el- 
ements are themselves (N, N) matrices. Similarly, we denote 
by ai an A^-element subvector of a. 

We now have 

ai = G U ix v (19) 
V 

where it must be understood that each term Cu'iv involves 
a matrix multiplication of order N. Matrix C has the sym- 
metry property Cu* = {Gin) . Especially, the diagonal ele- 
ments C\x are symmetric. 

We apply Crout's algorithm to Eq. II19II . We follow the 
procedure presented in Press et al. (1992), keeping in mind 
that instead of scalar elements we are now operating with 
matrices. 

We want to decompose C into 

C lk = Ujt$i (20) 

3 

where L\j = for j > I. Here L is a lower triangular matrix, 
whose elements are again (N, N) matrices. Note that the 
transpose sign refers to the element Lkj, not L itself. 
We write 

C ik = LikLlk + ^LijLlj (l>k) (21) 

j<k 

Cu = LuLJi + 'y ^LijLjj. (22) 

From this we can solve the elements of L, 

Uk = [L^iCl-Y^LkjLl)? (l>k) (23) 

i<k 

Lu = [Cu-J2 l H^ 1/2 - ( 24 ) 

3<l 

For each I, we first use Eq. 123H to solve Ljfc for k = 1 ... I — 1 
and then Eq. 1241 to solve La. Once L is known, elements 
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Xh can be solved by backsubstitution as 



Lkljak-y^LkjZj) k = l...L 
L~kk(zk -^2Lj k xj) k = L...l. 



(25) 
(26) 



3>k 



The procedure presented above contains operations be- 
tween (N, N) matrices and vectors of length N. These 
(N, N) matrices are nearly circulant, except for the corners, 
where they do not 'wrap around' like circulant matrices do. 
However, the deviation is small, and we may treat the ma- 
trices as circulant. 

Circulant matrix operations are most conveniently car- 
ried out in Fourier space. We pad each vector with zeros up 
to the next power of two, and use the Fast Fourier Transform 
(FFT) technique to perform the Fourier transforms. 

To each elementary operation involving a circulant ma- 
trix corresponds an operation in Fourier space. The corre- 
sponding Fourier operations are the following. 

(i) Matrix multiplication - element by element multipli- 
cation of Fourier transforms 

(ii) Matrix transpose - complex conjugate of the Fourier 
transform 

(iii) Square root of a matrix - Square root of the Fourier 
transform. 

(iv) Matrix inversion - inverse of the Fourier transform. 

The procedure on determining x = a can be sum- 
marized as follows. First perform Fourier transform to the 
circulant ring index i — i' of C a . Then carry out Cholesky 
decomposition in index I as described above, and store the 
resulting L matrix. 

For each vector a, carry out a Fourier transform along 
the ring index i, perform backsubstitution as given by 12511 - 
12611 . and do an inverse Fourier transform to obtain x. 

The total operation count of the above procedure is 
proportional to L 3 N\nN, as contrasted to L 3 N 3 of normal 
matrix inversion. Furthermore, the decomposition can be 
done 'in place' in the space of L 2 N elements, instead of 
L 2 N 2 . 



3.3 Covariance of the output map 

The covariance of the output map of MADAM is given by 



[P T (C„ + FC a F T )- 1 P]-\ 



(27) 



assuming that the noise is well modelled by the noise model 
The inverse of 1271 can be put into the form 

C" 1 = (28) 
P T C- 1 P - P T C- 1 F(C- 1 + F T C- 1 F)- 1 F T C- 1 P. 

We presented in Section [3.21 a procedure for evaluating 
C^a for arbitrary a. By running this procedure L times 
one can compute the inverse of matrix C a . Matrix C„ 1 + 
F T C~ 1 F can then again be decomposed and stored using 
the same procedure. When this is done, one can then easily 
compute any element of C^ 1 using formula 112811 . 



4 COVARIANCE OF THE COMPONENT 
AMPLITUDES 

4.1 General 

In this section we consider the computation of the covariance 
matrix C a . 

First we define reference values for the amplitude vector 
a. Suppose the actual coadded noise TOD, denoted by u, is 
known. We consider here the correlated noise component 
only. We fit the noise model Fa to the noise stream. A least- 
squares fit gives 



a= (F T F) _1 F T u. 



(29) 



Eq. 1291 gives the reference values or best estimates for the 
amplitude vector. We compute the covariance matrix as 



C a = (aa T ). 



(30) 



Let now y p be the original, uncoadded noise stream. We 
assume that the noise is stationary and its auto-correlation 
Cp-p' = {y P y P >} is known. 

The coadded noise stream is 



M-l 

-y 

m = 



(31) 



Here n is the number of samples per ring and M is the 
number of coadded circles (M = 60 for the nominal Planck 
scanning strategy). Index i = • • • N — 1 labels rings, j = 
• • • n — 1 labels samples on a ring, and m = • • • M — 1 
labels circles coadded into a ring. 

Let Fij be the values of the base functions I — 1 • • • L 
on a ring. We assume orthogonality and write 

p v = ( 32 ) 

The reference values for component amplitudes are, accord- 
ing to 1291 1 . given by 



an 



n—l 

E 

3=0 



n— 1 M-l 



^— ' M 

j=0 m=0 



(33) 



For uniform baselines, for instance, the reference value is 
simply the average of the noise over the ring. 

Next we calculate the theoretical covariance of the ref- 
erence values 13311 . The elements of the covariance matrix 
are given by 



Ca,ii'W — (5ij5i'(') 



(34) 



n— 1 

£ 

3,j'=0 
n—l 



PljPl 'i'Tf2 £ (y[inM+mn+j]V[i'nM+m'n+j']) 



m,m f =0 



M-l 



- 2^1 l * l 'j'~M2 C [{i-i')nM + {m-m')n+3-3'}- 

j,j'=Q m,m / =0 

The sum over m, m! can be combined into one sum over 
m" = m— m! . 



n—l M 
3,j'=0 m"=—l 



(35) 



M — \m"\ 



~C[(i — i , )nM-\-m"n J rj—j / ] • 
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This is a general formula for the elements of the covariance 
matrix C a . 



4.2 Exponential expansion of the auto-correlation 
function 

Covariance l|35[l is rather heavy to evaluate computationally, 
since it contains a three-dimensional sum. However, if the 
auto-correlation function can be expanded in exponential 
functions, the covariance can be computed in a very efficient 
way. This holds, e.g., for 1// and l// 2 noise. 

Suppose the auto-correlation function can be expanded 



c(t) = ^2b k c k (t) 

k 

where 

Cfc(i) = exp(-g k \t\) 



(36) 



(37) 



where gk is a selected set of characteristic frequencies and 
and bk are coefficients to be determined. 

We denote by C fc the covariance matrix that corre- 
sponds to an exponential auto-correlation function of the 
form I37H . Once the component covariances C h and coeffi- 
cients bk have been determined, the total covariance C a can 
be computed as 



2>c\ 



(38) 



The component covariances C k can be computed in a very 
efficient way making use of the basic properties of the expo- 
nent function. 

Expanding the auto-correlation as I36H is equivalent to 
expanding the power spectrum as 



(39) 



The expansion does not exist for arbitrary noise spectra. 
In Appendix A we calculate the expansion explicitly for a 
power-law spectrum of the form 



/>(/) = ( ?f) j 



(/ > /mm). 



(40) 



for — 2 ^ 7 < 0. Here f^ n is the knee frequency, at which 
the spectral power equals the white noise power a 2 / f s , a is 
the white noise std, and / s is the sampling frequency. 

For 1// noise (7 = —1) the expansion is particularly 
simple. If the characteristic frequencies gk are chosen uni- 
formly in logarithmic scale, the correct spectrum is obtained 
with 



fa 



(41) 



where A is the logarithmic interval in g k - 

The 1/f spectrum, as given by the expansion 11390 with 
coefficients 1411 . is shown in Fig. The 1/f form holds 
inside the frequency range / m j n < / < /max spanned by the 
characteristic frequencies gk- Below / m in the spectrum levels 
out, as can be seen from the figure. 

Another simple example is the l// 2 spectrum (7 = —2). 
In that case the desired spectrum is given by one single g 
component with coefficient 




Figure 1. 1/f noise power spectrum, as given by expansion 1391 
and Eq. 1411 . The solid line shows the pure 1/f noise spectrum. 
The dashed line shows the white noise level (P wn = c 2 / fa)- The 
dash-dotted line shows the total noise spectrum, including both 
components. The parameters were cr = 2700 fiK, / s = 76.8 Hz, 
/ kn = 0.1 Hz, A = 1, / min = 10" 5 Hz, / max = 10 Hz. 



2^ 2 / k 2 



(42) 



9 U 
and g = 27r/ mill . 

4.3 Component covariance matrices 

In this section we give explicit formulae for the elements 
of the component covariance matrix C fe that corresponds 
to the auto-correlation function 13711 . Derivation is given in 
Appendix A. Here we just quote the results. 

We use again the index notation, where indices i, 1 re- 
fer to rings and 1,1' to base functions. The elements of the 
component covariance matrix are given by 



Cunr = GkS£ t S M , exp ( - -^(i - i 



Cu'w = G k S kl S+, exp ( - - 



1)M) 
i - 1)M) 



(1 > 1 



(1 < 1 



(43) 
(44) 

(45) 



+G k(S+S-,+S k - l S+,). (i = i') 



Here / sp = f s /n, where n is the number of samples on a 
ring, represents the spin frequency of the instrument. If no 
coadding is applied, n can be chosen freely, and / sp does not 
need to have any connection to the scanning pattern of the 
instrument. In that case / sp represents simply the inverse of 
the chosen baseline length. 

Factors S + and S~ are defined as 



S+=VF y exp(-fi) 
— : Js 



3=0 
-1 



S k i =y^i ? !jexp(-^(n- 



-]))■ 



(46) 



(47) 
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Figure 2. Covariance between uniform baselines, as a function 
of the distance between rings. The solid line shows the theoretical 
curve, computed using Eqs. 1431 - 1451 , Diamonds present the val- 
ues determined from simulated noise. The noise parameters were 
a = 2700 ,uK, /„ = 76.8 Hz, / kn = 0.1 Hz, A = 1, / min = 10~ 5 
Hz, / max = 10 Hz. 



Gk = 



and 



Coadding brings in the factors 
1 [1 - exp(- 



M 2 [1 - exp( 



-f-M)f 

9k \]2 ' 

fsp >> 



Gl = 1 



M 1 



exp(-^) 



1 - 



1 1 -exp(-j±M) 



M 1 - exp(- 



9k ' 
fsp' 



(48) 



(49) 



If no coadding is done (M = 1), we have G k = 1 and G° k = 0. 
Factor Aij is defined by 



j-i 

Ay = VF^expC-f (j 



j'=o 



(50) 



and can be computed rapidly using the recurrence relation 
Ay = +F,j_i)exp(-^), (51) 



with the starting value A;o = 0. 

Formulae 1431 - ll45| l are very fast to evaluate numeri- 
cally, as compared to the general formula 135> . 

Figure [5] presents the theoretical covariance, computed 
using expansion 13811 . between uniform baselines for 1// 
noise with /k n = 0.1 Hz. Other parameters used were 
f s = 76.8 Hz, / gp = 1/60 Hz, n = 4608 and M = 60. 
We show in the same figure the covariance as determined 
from simulated 1// noise. We generated 10 realizations of 
noise TOD of one year length, and computed their auto- 
correlation using the Fourier technique. The agreement is 
very good. 

As another example we show in Table [5] the first ele- 
ments of the covariance matrix for Fourier components. In- 
dex / = 1 refers to the uniform baseline and indices I — 2 
and I — 3 (I — 4 and I = 5) to the sine and cosine of the 
first (second) Fourier mode, respectively. We have normal- 
ized all components to ^ . Jy = n. Elements of matrix F are 



Table 2. The first elements of the covariance matrix C a (in 
(^tK) 2 ). The noise parameters were the same as in Fig. [5] We used 
Fourier components as base functions, normalized as y\ F 2 = n. 
We show elements 1,1' = 1 ... 5 and i — i' = . . . 3. Index value 
/ = 1 refers to the uniform baseline and values I = 2 and I = 3 
(1 = 4 and I = 5) to the sine and cosine of the first (second) 
Fourier mode, respectively. The noise parameters were the same 
as in Fig. [5] 



/ 


i-i' 


V = 1 


I' = 2 


V = 3 


V = A 


V = 5 


1 





56049 





-2.41 





-0.668 




1 


31324 


-89.6 


1.10 


-44.9 


0.307 




2 


18707 


-28.9 


0.0531 


-14.5 


0.0133 




3 


12928 


-16.3 


0.0210 


-8.17 


5.25e-3 


2 








161 





1.64 







1 


89.6 


-1.42 


0.209 


-0.745 


0.0696 




2 


28.9 


-0.0751 


2.35e-4 


-0.0376 


5.88e-5 




3 


16.3 


-0.0297 


5.90e-5 


-0.0148 


1.48e-5 


3 





-2.41 





158 





-0.123 




1 


1.10 


-0.209 


0.133 


-0.139 


0.0617 




2 


0.0531 


-2.35e-4 


1.14e-6 


-1.18e-4 


2.85e-7 




3 


0.0210 


-5.90e-5 


1.73e-7 


-2.95e-5 


4.32e-8 


4 








1.64 





79.9 







1 


44.9 


-0.745 


0.139 


-0.401 


0.0524 




2 


14.5 


-0.0376 


1.18e-4 


-0.0188 


2.94e-5 




3 


8.17 


-0.0148 


2.95e-5 


-7.42e-3 


7.38e-6 


5 





-0.668 





-0.123 





79.0 




1 


0.307 


-0.0696 


0.0617 


-0.0524 


0.0334 




2 


0.0133 


-5.88e-5 


2.85e-7 


-2.94e-5 


7.13e-8 




3 


5.25e-3 


-1.48e-5 


4.32e-8 


-7.38e-6 


1.08e-8 



thus Fij = 1, F 2j = V2sm(2nj/n), F 3j = y/2 cos(2n j/n), 
Fij = \/2 sin(47rj'/ra), F$j = v2 cos(47rj/n), and F is given 
by Fij = Fij/n. We see that the dominant elements are those 
corresponding to uniform baselines. 



5 SIMULATION RESULTS 
5.1 Data sets 

We have produced two sets of simulated TOD. We refer to 
them as 'coadded' and 'uncoadded' data sets. 

The coadded data set mimicks the one year TOD from 
one Planck LFI 70 GHz detector. The scanning pattern was 
the following. The spin axis remained in the equatorial plane 
and was turned 2.4 arcmin every hour, so that after 8640 
hours the spin axis had turned 360 degrees. The detector 
turned around the spin axis with an opening angle of 85 
deg and spin frequency / sp = 1/60 Hz. The sampling fre- 
quency was f s = 76.8 Hz. We coadded data of 60 consecu- 
tive spin circles to form a ring. Our total data set consisted 
of 8640 rings, with 4608 samples in each. The sky coverage 
was 0.9964. 

The uncoadded data set was generated with a quite sim- 
ilar scanning pattern. The main difference was that instead 
of moving in steps, the spin axis turned continuously at the 
rate of 360 degrees in 8640 min. The sampling and spin fre- 
quencies as well as the opening angle were the same as in 
the first data set. Because the spin axis moved continuously, 
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consecutive circles did not fall on top of each other, and no 
coadding was done. The total length of the TOD was the 
same as in the coadded data set, i.e. 8640x4608 samples. 
The amount of data was equivalent to 6 days of one detec- 
tor Planck data, spread over the whole sky. This scanning 
pattern was rather artificial, but our purpose was only to 
demonstrate the use of MADAM in the case of uncoadded 
data. Full-scale simulations with realistic uncoadded one- 
year Planck data are beyond the scope of this paper. 

The underlying CMB map was created by the Synfast 
code of the HEALPix package (Gorski et al. 1999, 2004), 
starting from the CMB anisotropy angular power spectrum 
computed with the CMBFAST 3 code (Zeljak & Zaldarriaga 
1996) using the cosmological parameters fltot = 1.00, Oa = 
0.7, Q b h 2 = 0.02, h = 0.7, n = 1.00, and r = 0.0. We created 
the input map with HEALPix resolution iV s id e = 2048 and 
with a symmetric Gaussian beam with full width at half 
maximum (FWHM) of 14 arcmin. We then formed the signal 
TOD by picking temperatures from this map. Our output 
maps have resolution parameter A' s i < j c =512, corresponding 
to an angular resolution of 7 arcmin. 

We used the Stochastic Differential Equation (SDE) 
technique to create the instrument noise stream, which we 
added to the signal TOD. We generated noise with power 
spectrum 

p w = ( 1+ y)t' (/>/min) (52) 

with parameters a = 2700 fiK, knee frequency /k n = 0.1 Hz, 
and / min = 10~ 5 Hz. The white noise level 2700 fxK (CMB 
temperature scale) corresponds to the estimated white noise 
level of one 70 GHz LFI detector. We used the same noise 
spectrum for both data sets. 

We run our code on one processor of an IBM eServer 
Cluster 1600 supercomputer. 

5.2 Results for coadded data 

We show our results for the coadded data set in Tables 13151 
As a figure of merit we have used the rms of the residual 
noise map. The residual noise map was computed by sub- 
tracting from the output map a reference map. The reference 
map was created by coadding the pure signal TOD into a 
map. We then subtracted the monopole from the residual 
map and computed its rms value. 

The results for different numbers of base functions are 
given in Table [3] The given rms values are averages over 
10 noise realizations. We tried two sets of base functions: 
Fourier components and Legendre polynomials. Because a 
Fourier fit always includes an uniform baseline plus an equal 
number of sine and cosine functions, the total number L of 
base functions is always an odd number. The rms values 
continue improving when we increase the number of base 
functions. Fourier components give lower rms values than 
Legendre polynomials for the same number of components. 
In the rest of the simulations in this section we fitted only 
Fourier components. 

We show also the number of iterations and total CPU 
time taken by one run in the case of Fourier components. 

3 http://www.cmbfast.org 



Table 3. Average rms of the residual noise map for different 
numbers of base functions, for coadded data. We fit Fourier com- 
ponents and Legendre polynomials. The averages are taken over 
10 noise realizations. We show also the number of iteration steps 
and CPU time in the Fourier case. 





Legendre 


Fourier 






L 


rms/ fiK 


rms//iK 


iter 


CPU/s 


1 


110.634 


110.634 


16 


59 


2 


110.560 








3 


110.525 


110.485 


20 


65 


1 


110.481 








5 


110.451 


110.422 


28 


106 


7 


110.413 


110.386 


28 


12.3 


9 


110.387 


110.363 


32 


164 


11 


110.368 


110.347 


32 


191 


15 


110.344 


110.326 


36 


267 


25 


110.312 


110.301 


10 


530 


35 


110.298 


110.290 


11 


758 


45 


110.290 


110.283 


52 


1154 


65 




110.277 


56 


1959 







+ 


+ 
* 

L=3 


+ 


* 

+ : 

+ * : 
* 


* 

L=5 




* 

+ L=15 








- +"* 
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Figure 3. Rms of the residual noise map plotted against the 
inverse of the number of base functions, for Legendre polynomials 
(+) and Fourier components (*). The numerical values are given 
m Table El The std seems to converge towards value 110.26 /iK 
at the limit L — > n. 

Since the CPU time naturally depends on the computer used 
and may vary from run to run, the values quoted should not 
be taken too seriously. They merely give an idea how the run 
time increases with increasing number of base functions. 

In Fig. [3]we have plotted the rms values against the in- 
verse of the number of base functions. At the limit 1/L — > 
the rms values seem to be approaching the value 110.26 fiK. 
We expect that to be the std of the minimum-variance map 
(Section l2.2ll . The expected contribution from white noise to 
the residual map rms is 108.95 /j,K. This value was computed 
from the white noise sigma and the known distribution of 
measurements in the sky. 

We have compared results of fitting uniform baselines 
using MADAM and ordinary destriping technique. We got 
the destriping results by running MADAM with 1 = 0. 
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Table 4. Effect of misestimating the noise spectrum. We show the 
average residual map std (in fiK) for different numbers of Fourier 
components and for different assumed knee frequencies f^ n (upper 
panel), spectral slopes 7 (middle), and minimum frequencies / m i n 
(lower panel). The other two noise parameters were kept at their 
correct values (/ kn = 0.1 Hz, 7 = —1.0, and / m j n = 10~ 5 Hz). 
The lowest rms value on each row is denoted by an asterisk. The 
correct parameter value is shown in boldface. 



/kn 



L 


0.03 Hz 


0.05 Hz 


0.1 Hz 


0.2 Hz 


1 


110.634* 


110.634* 


110.634* 


110.634* 


3 


110.492 


110.484* 


110.485 


110.492 


5 


110.436 


110.420* 


110.422 


110.441 


7 


110.409 


110.387 


110.386* 


110.411 


9 


110.392 


110.366 


110.363* 


110.393 


11 


110.382 


110.353 


110.347* 


110.380 


15 


110.369 


110.335 


110.326* 


110.363 


25 


110.355 


110.319 


110.301* 


110.344 



7 



L 


-0.6 


-0.8 


-1.0 


-1.2 


1 


110.640 


110.635 


110.634* 


110.634 


3 


110.490 


110.483* 


110.485 


110.488 


5 


110.428 


110.419* 


110.422 


110.428 


7 


110.394 


110.384* 


110.386 


110.393 


9 


110.373 


110.361* 


110.363 


110.369 


11 


110.359 


110.346* 


110.347 


110.352 


15 


110.340 


110.326* 


110.326* 


110.330 


25 


110.321 


110.303 


110.301* 


110.304 


35 


110.314 


110.293 


110.290* 


110.292 






/min 






L 


10" 6 Hz 


10 5 Hz 


10~ 4 Hz 


10- 3 Hz 


1 


110.634* 


110.634* 


110.635 


110.655 


3 


110.485* 


110.485* 


110.485 


110.510 


5 


110.422* 


110.422* 


110.422 


110.463 


7 


110.386* 


110.386* 


110.386 


110.433 


9 


110.363* 


110.363* 


110.363 


110.415 


11 


110.347* 


110.347* 


110.347 


110.402 


15 


110.326* 


110.326* 


110.326 


110.386 


25 


110.301* 


110.301* 


110.302 


110.367 



At this limit the method reduces to pure destriping. The 
destriping result was 110.63444 (110.63443 fiK with 
MADAM). This indicates that when fitting uniform baselines 
only, the covariance plays very little role, but the baselines 
can be determined from the data alone with good accuracy. 

Keihanen et al. (2004) showed that fitting Fourier com- 
ponents beyond the uniform baseline with the ordinary de- 
striping technique, without using the covariance matrix, did 
not improve the results. In this work we have found a clear 
improvement. This indicates, that information about the 
noise spectrum is important when fitting base functions 
other than the uniform baseline. 

We have also studied the effect of misestimating the 
noise spectrum. We varied in turn each of the three noise 
parameters (knee frequency, spectral slope, and minimum 
frequency) while keeping the other two at their correct val- 
ues (/ kn = 0.1 Hz, 7 = -1.0, fmin = KT 5 Hz). We then 
recomputed the covariance matrix C a with the new param- 
eter values and rerun the map estimation. We fitted Fourier 
components only. The results are shown in Table 2] 



Table 5. Pixclization noise. Rms of the residual map for noise- 
free TOD. The error is due to finite pixelization of the sky. 



L 


rms/ /jK 


1 


0.138 


3 


0.207 


5 


0.233 


7 


0.257 


9 


0.274 


11 


0.289 


15 


0.319 


25 


0.359 


35 


0.380 


45 


0.393 



It is perhaps surprising that for small L underestimat- 
ing the knee frequency or assuming a less deep slope seems 
to improve the results. This can be understood as follows. 
When L is small, the noise is not perfectly modelled by the 
base functions. There is an error involved, related to the 
higher Fourier components that have not been included in 
the analysis. This error affects the estimation of the lower 
components, leading to a too high variation in their ampli- 
tudes. The error in the covariance matrix, caused by mises- 
timation of the noise spectrum, partly compensates for this 
error. We notice that the best results are obtained with a 
spectrum (less deep a slope or lower knee frequency) which 
leads to a smaller covariance for the low-frequency Fourier 
components. Smaller covariance tends to restrict the varia- 
tion of the amplitudes, thus decreasing their error also. With 
larger L the phenomenon disappears, and the lowest rms is 
obtained with the correct noise spectrum, as expected. 

Table shows results from a run with noise-free data. 
The TOD contained only the contribution from the CMB 
signal, but no instrument noise. The error that still remains 
in the map is due to 'pixelization noise' (Dore et al. 2001) 
caused by the finite size of sky pixels. The pixelization error 
increases with increasing number of base functions, but is 
very small compared with the error due to instrument noise. 

5.3 Results for uncoadded data 

If no coadding is applied, the length n of a ring is not de- 
termined by the scanning pattern of the instrument, but is 
a free parameter to be chosen at will. We then have two pa- 
rameters to select: the number L of base functions and their 
length n. 

To keep things simple, we tried two schemes. First we 
kept the baseline length fixed at one minute and varied the 
number of Fourier components that we fitted. Secondly, we 
fitted uniform baselines only (L = 1) but varied their length. 

Results from the first case are shown in Table |5] The 
baseline length was fixed at n — 4608 samples (one minute). 
We show again the average rms of the residual noise map, 
averaged over 10 realizations of noise. The white noise level 
is higher than in the coadded case by a factor of \/60- The 
expected white noise contribution to the map rms is 844.0 

(J.K. 

Table |7| shows results of fitting uniform baselines of dif- 
ferent lengths. The first column gives the length n of the 
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Table 6. Average residual noise rms for different numbers of 
Fourier components, for the uncoadded data set. The ring length 
was n = 4608. Also shown are the number of iteration steps and 
the CPU time for one run. 



L rms/^tK iter CPU/s 



1 


857.131 


28 


68 


3 


855.782 


36 


102 


5 


855.294 


18 


161 


7 


838.373 


52 


205 


9 


854.842 


56 


280 


11 


854.718 


60 


336 


15 


854.555 


64 


436 


25 


854.357 


76 


847 



Table 7. Results of fitting uniform baselines of different lengths 
to uncoadded data. The first three columns give the baseline 
length as the number of samples and in seconds, respectively, 
and the number of baselines per one minute of TOD. The next 
columns give the average residual noise rms, number of iteration 
steps and total CPU time taken by one run. 



n 


t/s 


N/min 


rms / /iK 


iter 


CPU/s 


4608 


60.0 


1 


857.131 


28 


68 


2304 


30.0 


2 


856.220 


32 


80 


1152 


15.0 


i 


855.837 


36 


96 


576 


7.5 


8 


855.202 


48 


123 


288 


3.75 


16 


854.769 


64 


150 


144 


1.88 


32 


854.463 


84 


215 


72 


0.94 


64 


854.271 


116 


331 


36 


0.47 


128 


854.169 


160 


672 


18 


0.23 


256 


854.116 


224 


1284 


9 


0.12 


512 


854.092 


320 


2889 



baseline, as the number of samples. The second column gives 
the baseline length in seconds. The third column shows the 
number of baselines per minute (4608/n). The shortest base- 
line we tried consisted of only 9 samples. 

The third column of Table Q an d the first column of 
Table |S| are comparable, since they give the total number of 
unknows per one minute of TOD. We see that for a given 
number of unknowns, fitting Fourier components works bet- 
ter than fitting uniform baselines. However, when we com- 
pare CPU times, we see that fitting uniform baselines is 
more effective. 

As in the case of coadded data, we compared results 
of fitting uniform baselines using MADAM and ordinary de- 
striping technique. The results are shown in Table |H] With 

Table 8. Comparison between MADAM and destriping. 









MADAM 


destriping 


n 


t/s 


N/min 


rms//iK 


rms / /iK 


4608 


60.0 


1 


857.131 


857.135 


2304 


30.0 


2 


856.220 


856.239 


1152 


15.0 


4 


855.837 


856.779 


576 


7.5 


8 


855.202 


861.118 


288 


3.75 


16 


854.769 


875.798 



one minute baselines the difference between the methods is 
small, but increases with decreasing n, in favour of MADAM. 
The rms value obtained with destriping reaches a minimum 
around 0.5 min baseline length, while with MADAM the val- 
ues continue improving. With small values of n MADAM is 
clearly superior to basic destriping. 



6 CONCLUSIONS 

We have presented a new map-making method for CMB 
experiments called MADAM. The method is based on the 
well known destriping technique, but unlike basic destriping, 
it also uses information on the known statistical properties 
of the instrument noise in the form of the covariance matrix 
of the base function amplitudes. We have shown that with 
this extra information the CMB map can be estimated with 
better accuracy than with pure destriping. 

We have tested the method with simulated coadded 
Planck -like data. As a figure of merit we have used the 
rms of the residual noise map. Our simulations show that 
fitting more base functions clearly improves the accuracy of 
the output map, with the cost of increasing requirements for 
CPU time and memory. 

We have shown theoretically that the map esti- 
mate given by MADAM approaches the optimal minimum- 
variance map when the number of fitted base functions in- 
creases. In practice it is not possible to reach the exact 
minimum-variance map using MADAM, due to CPU time 
and memory limitations. Still, MADAM provides a fast and 
efficient map-making method. By varying the number of 
base functions the user may flexibly move from a very fast 
but less accurate map-making (small L) to a more accurate 
but more time-consuming map-making (large L), depending 
on what is desired. 

We also demonstrated the use of MADAM for uncoad- 
ded data. Although the data set used was quite artificial, 
in the sense that it does not mimick data from any existing 
CMB experiment, the method was shown to work well for 
uncoadded data also. 

The current implementation of the method is a serial 
one. With a parallelized version we expect to be able to 
process data sets equivalent to full-year uncoadded Planck 
data. 
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The total power spectrum corresponding to the auto- 
correlation iAll is given by 



9l + (2^/) 2 ' 



(A5) 



We pick the frequencies gt uniformly in logarithmic 
scale inside some range [/min, /max] and use the ansatz 

b k = Ag2 +1 A (A6) 

where A = \n(gk+i/ gk) is the logarithmic step in g and A 
is a constant to be determined. The total power spectrum 
becomes 



9l + (2^/) 2 



(A7) 



We transform the sum into an integral 

f/max „ 7+ 2 , 



P(f) = A 



2(f ! 



/min 



g* + (2tt/)2 g 



(A8) 



r°° 2o 7+1 
~ A J g 2 + (27r/) 2 dg (/min < f < /max) ' 

The integral converges for — 2 < 7 < 0, 



i[( 7 + 2)tt/2] 



(A9) 



We choose 



APPENDIX A: COVARIANCE MATRIX FOR A 
POWER-LAW POWER SPECTRUM 

We discussed the computation of the covariance matrix C a 
in Section 4. In this appendix we present some technical 
calculations which were omitted there. 



Al Exponential expansion for a power-law 
spectrum 

Assume that the auto-correlation function of the noise can 
be expanded as 



c(t) = ^6fc exp(-g k \t\) 



(Al) 



where and gt is a selected set of characteristic frequencies. 
We now derive the coefficients bk in the case of a power-law 
power spectrum of the form 



p(f) oc r 



(A2) 



for -2 ^ 7 < 0. 

The power spectrum P(f) and the auto-correlation 
function c(t) of stationary noise are related by the cosine 
transform 



P(f) 



c(t)cos(2nft)dt. 



(A3) 



The auto-correlation function exp(— g\t\) corresponds to the 
power spectrum 

2ff 



P(f,g) 



9 2 + (2tt/) 2 



(A4) 



A = — -(2^/ kn )" 7 sin[( 7 + 2)n/2] 

Is 7T 

to obtain the desired power spectrum 

p(f) = C(i^r (/ min «/«/ m 

Js Jkn 



(A10) 



(All) 



Here a is the white noise std, / s is the sampling frequency, 
and /kn is the knee frequency, at which the power of the 
power-law noise component equals the white noise power 
cr 2 // s . The maximum / max should well exceed the knee fre- 
quency. 

The coefficients bk are given by 



b k = ^-(2^/ kn )- 7 sin[(7 + 2V/2]^ +1 A 

fa 7T 



(A12) 



for —2 < 7 < 0. Especially, for 1/f noise (7 = —1) we have 
the simple formula 



fs 



(A13) 



The case 7 = — 2 requires a separate treatment. We see 
directly from Eq. 1A4L that the desired spectrum is given 
by one single g component with coefficient 



6 _27r 2 / k y 2 

9 fs ' 

The spectrum then has the form 



P(f) = 



ft 



/min + f 2 fs 



(A14) 



(A15) 



where / m j n = g/(2w). The spectrum behaves as oc / 2 at 
/ 3> /min and levels out below /mm- 
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A2 Component covariance matrices 

Once the expansion has been found, the covariance 

matrix can be computed as 



(A16) 



In the following we calculate the component covariance ma- 
trices C fe . 



In Section 4 we derived the general formula 



(A17) 



n-l M-l 

E Ftifi'j'JtfE E C[(i-i>M+(m-m')n+j-j']' 
j,j' = m,m'— 

Assume then that the auto-correlation function is of the 
exponential form 

(Mp') = Ck, P - v > = exp(-g k \t\) = cxp ( - (A18) 

where / s is the sampling frequency and indices p,p' label 
samples along the TOD. We substitute this into Eq. IA17II . 
The covariance becomes 



n-l M-l 

Cu'W = 51 ^^'Vjjfj E 

jj'=0 m,m' — 



(A19) 



x 



exp ( — ^r-|(i — i')nM + (m — m')n + j — 



We treat cases |i — i | > and i = i separately. 

(i) \i — i'\ > 0. If i — i' > the quantity inside the brackets 
is positive, and we can split the four-dimensional sum into 
a product of four sums as 

C%nv = exp(-f (i-i'-l)nM) 

M-l M-l 

X M?E e 77 E e 77 

m— m=0 



re-1 



n-l 



x ]TiV~^' E f iV e-^("- J '>. (A20) 

3=0 j'=0 

We have arranged the terms in such a way that the argument 
of an exponent function is always negative. This is helpful 
in numerical evaluation. The sum over m, m! can be carried 
out analytically, yielding 

M-l M-l 
1 \ ~* -Sk mn \ -2k(M-l-m)n 

Gk = a^E 6 /a E e /B 

m — m=0 



1 (1-e «p ) 
(1 — e /s p) 2 



(A21) 



where / sp = f s /n. 

The elements for which i — i' < are obtained from the 
symmetry relation,^,,,, = C$ m . 

(ii) The case i = i' is more complicated, since the quantity 
inside the brackets in Eq. <A19t takes both positive and 
negative values. We split the sum into three terms (m = m' , 
m > m' , and m < m') and the m = m' term further into 



three terms (j = f , j > f , and j < f) 



n.— l m-l 

E Wwjji E e-^—'W> 

M—l 

1 -Si.{m-m' -l)n 

M 2 2-^i e 

m— m' <m 

7i — 1 n — 1 

x^Fye-^^F !V e-*^ 



3=0 
M-l 



3*'=0 



71 E E e-*(»'— ^ 



M 2 



m'— m<m' 
n-l 



))-l 



x E F ^ e ; 

3=0 3'=0 
n— 1 n— 1 



+ tjE'-'- • itE'vE'— * 



0-3') 



3=0 

n-l 



j'=0 j'<j 

L (i'-3) 



3'=0 3<j' 

The sum over m,m' can again be calculated analytically, 



G° k 



ivi — 1 

AJ2 E E exp(-g-(m-m'-l)) (A23) 



M 2 

m— m' <m 

1 1 

M l-exp(-|M 



1 - 



I l-exp(-J^M) 



M l-exp(-|M 

Formula EH may seem complicated, but is easy and fast 
to evaluate numerically. 

Equations HA20t and 1A22L when written in compact 
form, give the formulae 1431 - 145(1 in Section 4. 



